Chemical Reaction between Single Hydrogen Atom and Graphene 
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We study chemical reaction between a single hydrogen atom and a graphene, which is the 
elemental reaction between hydrogen and graphitic carbon materials. In the present work, classical 
molecular dynamics simulation is used with modified Brenner's empirical bond order potential. 
The three reactions, that is, absorption reaction, reflection reaction and penetration reaction, are 
observed in our simulation. Reaction rates depend on the incident energy of the hydrogen atom 
and the graphene temperature. The dependence can be explained by the following mechanisms: 

(1) The hydrogen atom receives repulsive force by 7r-electrons in addition to nuclear repulsion. 

(2) Absorbing the hydrogen atom, the graphene transforms its structure to the "overhang" 
configuration such as sp 3 state. (3) The hexagonal hole of the graphene is expanded during the 
penetration of the hydrogen atom. 
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I. INTRODUCTION 

Chemical reaction between hydrogen and graphite 
is a basic process of the chemical sputtering between 
plasma and a divertor plate in plasma confinement 



experiments 



1.2.3.4.5 



The divertor plate composed of 



graphite tiles or carbon fiber composites is bombarded 
with hydrogen plasma. Eroding the carbon wall, hy- 
drogen ions yield H2 and several sorts of hydrocarbon 
molecules, i.e., CH X ,C2H X and so on. The hydrocar- 
bon molecules misbehave as impurities for plasma con- 
finement experiments. In order to reduce the hydrocar- 
bon impurities in the plasma, we need to understand 
the erosion mechanism of carbon walls and the creation 
mechanism of hydrocarbon molecules. However, these 
mechanisms are not clarified yet. We, therefore, reveal 
the chemical reaction between hydrogen and graphite by 
computer simulation. 

In the present work, we study the chemical reaction 
between a single hydrogen atom and a graphene with 
classical molecular dynamics (CMD) simulation. It is 
reasonable to consider that ions of the hydrogen plasma 
combine with electrons and become neutral atoms be- 
fore interacting with the carbon wall. Therefore, we se- 
lect a neutral hydrogen atom as an injected particle. A 
graphene is the elemental component of graphitic car- 
bon materials.— The chemical reaction between the single 
hydrogen atom and the graphene, therefore, is regarded 
as the elemental reaction between hydrogen and various 
graphitic carbon materials. 

We measured only an absorption rate in our previous 
workji where multi-layer graphite was treated. In the 
present work, we evaluate a reflection rate and a pene- 
tration rate in addition to the absorption rate. We also 
obtain the incident hydrogen energy dependence and the 



graphene temperature dependence of each reaction rate. 
In <Jnl the simulation model and method are described. 
Simulation results are shown in ^IIIl We discuss the fea- 
ture of the chemical reaction between the hydrogen atom 
and the graphene in §IVI This paper is concluded with 
summary in j fvl In addition, we note the modification 
of Brenner's reactive empirical bond order potential in 
Appendix. 



II. SIMULATION METHOD 

In the present work, we adapt CMD simulation with 
the NVE condition, in which the number of particle, vol- 
ume and total energy are conserved. The second order 
symplectic integration^ is used to solve the time evolution 
of the equation of motion. The time step is 5 x 1CP 18 s. 
We use a modified Brenner's reactive empirical bond or- 
der (REBO) potential^ 
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where nj is the distance between the i-th and the j-th 
atoms. The functions V^-i and represent repulsion 

and attraction, respectively. The function bij generates 
multi-body force. We show details of the modified Bren- 
ner's REBO potential in Appendix. 

Figure [1] shows the present simulation model. The hy- 
drogen atom is injected into the graphene composed of 
160 carbon atoms. The center of mass of the graphene 
is set to the origin of coordinates. The surface of the 
graphene is parallel to the x-y plane. The size of the 
graphene is 2.13 nm x 1.97 nm. The graphene has no lat- 
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FIG. 1: Simulation model. There are 160 carbon atoms and 
an injected hydrogen atom. 



tice defects and no crystal edges due to periodic bound- 
ary condition toward x and y directions. The graphcnc 
temperature is defined by 



~ carbon 2 

E 



3Nk h 



2m, 



(2) 



where pi and Wj are the momentum and the mass of 
the i-th carbon atom, respectively. N is the number of 
carbon atoms and kb is the Boltzmann constant. The 
symbol ^^ arbon denotes the summation over the carbon 
atoms. The carbon atoms obey the Maxwell-Boltzmann 
distribution in the initial state of the simulation. 

The hydrogen atom is injected parallel to the z axis 
from z — 4 A. We repeat 200 simulations where the x 
and y coordinates of injection points are set at random. 
As a result, we obtain the histograms, which give reac- 
tion rates. The incident energy E\ determines the initial 
momentum pn(0) = (0, 0,prj) of the hydrogen atom as 



po = \j2mnEi, (3) 
where tor is the mass of the hydrogen atom. 



(a) 



(b) 




(c) 



(d) 




FIG. 2: Snapshots of the absorption reaction, (a) The hy- 
drogen atom with the incident energy of 3 eV is injected, (b) 
The hydrogen atom makes a covalent bond with the nearest 
carbon atom, (c) The nearest carbon atom is pulled out of 
the graphene and the "overhang" configuration appears, (d) 
The atoms are relaxed to the stable structure with oscillation. 



The nearest carbon atom is pulled out of the surface of 
the graphene as a sp 3 configuration (Fig. We call this 
phenomenon "overhang" . The hydrogen atom remains 
above the nearest carbon atom while oscillating. In the 
reflection reaction, the graphcnc reflects the incident hy- 
drogen atom to the region of z > 0. In the penetration 
reaction, the incident hydrogen atom passes through the 
graphene and goes away to the region of z < 0. It is 
observed that the graphene expands the hexagonal hole 
while the hydrogen atom is penetrating. 



III. RESULTS 



B. Incident energy dependence of reaction rates 



Three kinds of reactions between the single hydrogen 
atom and the graphene are observed in our CMD simu- 
lation. They are absorption reaction, reflection reaction 
and penetration reaction. The properties of the reactions 
are described in the following. 



A. Dynamics of three reactions 

In the absorption reaction, the hydrogen atom and the 
nearest carbon atom are bound by a new covalent bond. 



Figure [3] shows the incident energy dependence of each 
reaction rate in the case that the initial graphcnc tem- 
perature is 300 K. Three kinds of reactions dominate in 
different incident energy E\ respectively In the case of 
Ej < 1 eV, almost all of the incident hydrogen atoms are 
reflected. For 1 eV < E\ < 7 eV, the absorption reac- 
tion becomes dominant. The reflection reaction becomes 
dominant again for 7 eV < Ej < 30 eV. The penetration 
reaction behaves as the dominant process for E\ > 30 eV. 
It is also observed that the absorption rate has the small 
peak around 24 eV in Fig. [3J 
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FIG. 3: Incident energy dependence of the absorption, the 
reflection and the penetration rates. Dash-dotted line with 
open circle, long-dashed line with filled triangle, and short- 
dashed line with square denote the absorption, the reflection 
and the penetration rates, respectively. 



7 eV < Ei < 30 cV in Fig. H From this fact, it is 
deduced that two kinds of repulsive force work between 
the incident hydrogen atom and the graphene. To prove 
it, we plot the potential energy between the hydrogen 
atom and the graphene in Fig. [5) where the hydro- 
gen atom is located just above the nearest carbon atom 
at the distance w and the other carbon atoms are re- 
laxed. From Fig. we can confirm the existence of 
two kinds of repulsive force between the incident hydro- 
gen atom and the graphene. The first repulsive force for 
w < 1.0 A is due to the repulsive term in Eq. ([T]) 
and corresponds to nuclear repulsion. The second repul- 
sive force for 1.6 A < w < 1.8 A is derived from the 
multi-body force in the term bij in Eq. (TTJ. The ex- 
istence of the second repulsive force was also confirmed 
by ab-initio calculations . 1 1 1 1 2 It was considered that it— 
electrons over the graphene generate the second repul- 
sive force. The energy height of the potential wall of the 
second repulsive force is estimated to be about 0.5 cV. 
The hydrogen atom with the incident energy of 0.5 eV or 
more, therefore, can enter the region that w < 1.6 A, in 
which the other mechanism derives the absorption reac- 
tion and the penetration reaction. The mechanisms for 
E\ > 0.5 cV and w < 1.6 A are described in the sub- 
sequent subsections. Thus, the reflection rate starts to 
decrease at E\ ~ 0.5 eV in Fig. [3J 



C. Graphene temperature dependence of reaction 
rates 

We also investigate the initial graphene temperature 
dependence of reaction rates (Fig. rj}. As the initial 
graphene temperature rises, the absorption rate tends to 
broaden to the region of low incident energy. In con- 
trast, the reflection rate drops. However, the graphene 
temperature hardly affects the penetration rate. 



IV. DISCUSSIONS 

The three reactions were observed in the present simu- 
lation. The incident energy dependence and the graphene 
temperature dependence of the three reaction rates are 
also obtained. The similar incident energy dependence of 
three reactions was observed in the system of multi-layer 
graphite and a large amount of hydrogen atoms^ It is, 
therefore, important to understand the mechanism of the 
chemical reaction between a single hydrogen atom and a 
single graphene to argue atomic-scale processes in vari- 
ous systems composed of hydrogen and graphitic carbon 
materials. 



A. Two kinds of repulsive force 

It was observed that the reflection reaction dominates 
in two ranges of the incident energy, i.e., E\ < 1 eV and 



B. Absorption mechanism 

Figure [6] shows the potential energy contour in the 
(it, w) parameter space. Here it is the height of the near- 
est carbon atom from the surface of the graphene and w is 
the distance between the hydrogen atom and the nearest 
carbon atom as shown in Fig. [71 There is the minimum 
potential point at (u,w) = (0.5 A, 1.1 A), which corre- 
sponds to the "overhang" configuration. This analysis by 
the potential energy contour claims that the "overhang" 
configuration is the most stable state. 

The trajectory of the absorption reaction is represented 
by arrow ® in Fig. [5J The initial state corresponds to the 
point of (it, if) = (0 A, 4 A). The hydrogen atom and the 
graphene start interaction around (u,w) ~ (0 A, 1.8 A). 
The arrow ® shows that, with tumbling down the slope 
of the potential energy contour, the state falls into the 
potential minimum point (u,w) ~ (0.5 A, 1.1 A), which 
indicates the "overhang" configuration. 

Until the hydrogen atom overcomes the second repul- 
sive force, the graphene cannot transform its structure 
to the "overhang" configuration. Thereby, the incident 
energy of 0.5 eV, which is corresponding to the energy 
height of the potential wall of the second repulsive force, 
is the lower limit to occur the absorption reaction. The 
absorption rate rises from Ej ~ 0.5 cV in contrast to 
the reflection rate by the second repulsive force and has 
a peak at E\ = 3 eV. 



4 




FIG. 4: Graphene temperature dependence of the absorption, the reflection and the penetration rates. We varied the graphene 
temperature such as K, 300 K, 600 K, 800 K, 1000 K, 1500 K and 2000 K. 



> 

>, 
CD 

i 

o 
c 

CD 



c 

CD 
O 
0_ 




'.4 0.6 0.8 1 1.2 1.4 o 1.6 1. 
w : CH distance [A] 



FIG. 5: The potential energy between an incident hydrogen 
atom and a graphene. The potential energy is calculated by 
the modified REBO potential model. The position of the 
incident hydrogen atom is set to be (xo,yo,zo + w), where 
(xo,Vo,zo) is the position of the nearest carbon atom and 
w > 0. The other carbons are relaxed to the total potential 
energy minimum state for each w. 
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FIG. 6: Potential energy contour in the u and w parameter 
space. The parameters u and w axe indicated in Fig. [7] 
The potential energy is higher than 5 eV in the white area. 
The trident arrow represents the trajectory of each reaction. 
The numbers 1, 2 and 3 correspond to the absorption, the 
reflection, and the penetration reactions, respectively. 
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FIG. 7: Schematic picture of the "overhang" configuration. 
White and gray circles represent the incident hydrogen atom 
and the carbon atoms, respectively. Single lines represent 
covalent bonds. The parameter w is the distance between 
the incident hydrogen atom and the nearest carbon atom. 
The length u of the double line represents the height of the 
overhang from the surface of the graphene. 



C. Reflection mechanism 

In the reflection reaction for 7 eV < E\ < 30 eV, the 
incident hydrogen atom bounds back from the potential 
wall Vr^j in Eq. ([1]), which is represented by the region 

of w < 1.0 A in Fig. [5] and the white region in Fig. [6] 
After bounding, the hydrogen atom goes away from the 
graphene without connecting the nearest carbon atom. 
Therefore, the graphene keeps the flat sheet configuration 
and does not transform its structure to the "overhang" 
configuration. The trajectory of the reflection reaction 
by Vj^j is drawn as arrow (2) in Fig. [6] 

Here, to make clear a distinction between the absorp- 
tion reaction and the reflection reaction for Ej > 1 eV, 
we introduce two typical time lengths At u and At w . The 
time length At u is defined as the time length necessary 
for the graphene absorbing the hydrogen atom to trans- 
form its structure from the "flat sheet" configuration to 
the "overhang" configuration. Strictly speaking, At u de- 
pends on a lot of parameters, for example, the incident 
energy and the incident position of the hydrogen atom, 
the graphene temperature, and so on. But, to estimate 
the typical time length At u , we consider the only sim- 
ple overhang process, which is represented by the fol- 
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lowing trajectory in the parameter space Fig. [6j the 
configuration of the atoms is transformed from the start 
point that (u,w) = (0 A, 1.1 A) to the end point that 
(u, w) = (0.5 A, 1.1 A) along a straight line w = 1.1 A. 
We, moreover, assume that the initial velocities of all the 
atoms are zero. The potential function U a h along the 
above path is approximated by the following harmonic 
oscillator: 



(m c + m H )w 2 (w - u ) 5 



U , 



(4) 



where we use, as the mass, the sum of toc = 12.0 amu 
and mjj = 1-00 amu, because the hydrogen atom and 
the nearest carbon atom move as a rigid body in our 
assumption that if is fixed to the constant length of 
1.1 A. From Fig. [6l we have the potential minimum point 
Mo = 0.5 A, the minimum potential-energy Uo = 4.84 eV, 
and oj = 1.45 x 10 14 s _1 . Thus, we can estimate At u as 
follows: 



At u = t f — ) = 1-08 x 10" 1 
4 V uj 



(5) 



This estimated value of At u is comparable to the CMD 
simulation result At^ im - ~ 1.28 x 10~ 14 s, which was ob- 
tained under the condition that the degrees of the free- 
dom except the parameter u are fixed to the initial values. 

The other time length At w is defined as the time length 
in which the hydrogen atom can stay in the region that 
w < 1.6 A. To estimate At w , we adapt the alternative as- 
sumption that the hydrogen atom moves as a free particle 
for w > 0.9 A along the straight line of u — A and 
collides with the potential wall VSi at (u, w) 

= (o A, 

0.5 A). From this assumption and Eq. ([3]), At w is given 
by 



At, 



2/toh 
Po 



2mn 



(0) 



where / = (1.6 - 0.9) A = 0.7 A. From Eq. ©, it is 
obtained that Atu, is proportional to \/\fE[. On the 
other hand, Eq. © shows that At u is independent of 
Ei. 

Comparing these time length, we consider the follow- 
ing two cases. In the first case that At w > At u , the hy- 
drogen atom connects with the nearest carbon and the 
graphene transforms its structure to the "overhang" con- 
figuration, before the hydrogen atom escapes to the re- 
gion that w > 1.6 A. The hydrogen atom, therefore, is 
absorbed by the graphene. As the incident energy in- 
creases, At„, becomes smaller than At u . In this case 



At n > Atu, = l\/2mn/EiJ , the hydrogen atom escapes 

before the graphene traps the hydrogen atom. This pro- 
cess is regarded as the reflection reaction. We can derive 
the following condition necessary for the reflection reac- 
tion: 



Ei > 



2/ 2 m H 
At„ 2 



0.84 eV. 



(7) 



The incident energy which satisfies the condition Atu, = 
At u is estimated as Ef lm ~ 1 eV in our CMD simulation 
where the degrees of freedom except w are fixed to the 
initial values. By comparison between Ef lm and the con- 
dition Eq. (J7J) , it is considered that the our assumption 
is proper. In the above discussion, the hydrogen atom 
is located on the vertical axis over the nearest carbon 
atom. However, In the present simulation, the hydrogen 
atom seldom exists just above the nearest carbon atom, 
because the x and y coordinates of the incident hydro- 
gen atom are set at random. Thereby, the repulsive force 
by Vj^j becomes weaker than that of the potential en- 
ergy contour in Fig. [6l The time length At„, becomes, 
therefore, longer than the estimated value of Eq. j6|. 
The hydrogen atom which deviates from the vertical axis 
over the nearest carbon atom needs higher incident en- 
ergy than the estimated value of Eq. (J7J. Consequently, 
the incident energy of 0.84 eV in Eq. (J7J is the lower 
limit to occur the reflection reaction by . 



D. Penetration mechanism 

We describe the dynamics of the penetration reaction. 
We notice for the present simulation that the graphene 
expands the hexagonal hole during the penetration of 
the hydrogen atom. Figure [5] shows the potential en- 
ergy contour with two parameters, i.e., the distance w 
and the length v of the side of the hexagonal hole (See 
Fig. [9|). We note that the hydrogen atom is located 
above the center of the hexagonal hole unlike the layout 
of Fig. [Jj The C-C bond length of the stable graphene 
structure is 1.42 A. The interaction force acts on the hy- 
drogen atom and the graphene in w < 1.11 A. There is 
the potential minimum region of eV in the area that 
v = 1.42 A and w > 1.11 A, which is the incident state 
of the hydrogen atom. If the size of the hexagonal hole v 
is fixed to 1.42 A, the energy height of potential wall is 38 
eV at (v,w) = (1.42 A, A). In this case, the hydrogen 
atom needs the incident energy of 38 eV or more to pen- 
etrate the graphene. However, the penetration reaction 
with the incident energy of less than 38 eV is observed 
in the present simulation result Fig. [3l The difference 
between the estimation and the simulation result is ex- 
plained by the expansion mechanism of the hexagonal 
hole of the graphene. If carbon atoms move along the 
bottom of the potential energy valley in Fig. [8j the pa- 
rameter v increases from 1.42 A to 1.58 A with decreasing 
w. Thus, the hexagonal hole is expanded as the hydro- 
gen atom approaches the graphene. As a consequence, 
the energy height of the potential wall is lowered to 13 
eV at (v, w) = (1.58 A, A). This expansion lets the hy- 
drogen atom penetrate in the incident energy of less than 
38 eV. 

Here, we indicate that the carbon atoms can expand 
the hexagonal hole before reflecting the hydrogen atom 
with the incident energy of 13 eV. We define At' w as 
the time length for the hydrogen atom to approach the 
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graphene from w = 1.11 A. The time length At' w is given 
by 



AC = ^=I\/2*=2.18X10-»B. (8) 
Po V 2E i 

where I' = 1.11 A and Ej is set to 13 eV. On the other 
hand, the potential energy around v = 1.58 A, where w is 
fixed to A, is approximated by the following harmonic 
oscillator: 



^holc(w) 



(9) 
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where mhoic = 6toc, = 1-58 A. Uq = 13 eV and to' = 
5.21 x 10 14 s^ 1 from the potential energy contour Fig. 
[8l For this approximation, we obtain the time length 
At v to accomplish the expansion of the hexagonal hole 
as follows: 

At v = i f = 3.01 x 10~ 15 s. (10) 

Both At' w and At v are on the same order of femtosecond. 
In addition, the time length At' w of Eq. (|SJ) becomes prac- 
tically larger than 2.18 x 10~ 15 s because of deceleration 
due to repulsion. Therefore, the carbon atoms can ex- 
pand the hexagonal hole in response to the approach of 
the hydrogen atom. 

Next, we consider the small peak of the absorption rate 
at Ei = 24 eV, at which the hydrogen atom has enough 
incident energy to penetrate the graphenc. The incident 
energy of the hydrogen atom diffuses into the graphenc. 
Therefore, the hydrogen atom has no longer the necessary 
incident energy to escape from the graphene. From the 
energy diffusion, it is understood that the peak of the 
absorption reaction at 24 eV is caused by the hydrogen 
atom absorption on the reverse side of the graphenc. The 
absorption reaction on the reverse side was confirmed in 
the present simulation. As long as the hydrogen atom 
is absorbed, the graphene transforms its structure to the 
"overhang" configuration where the nearest carbon atom 
is pulled into the reverse side of the graphenc. 



E. Graphene temperature dependence of reaction 
rates 



FIG. 8: Potential energy contour in the v and w parameter 
space. The parameters v and w are indicated in Fig. [9] The 
Interaction force does not act on the hydrogen atom and the 
graphene when the state is located on w > 1.11 A. 
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FIG. 9: Expansion of the hexagonal hole of the graphene. 
White circle represents the incident hydrogen atom. Gray 
circles represent carbon atoms. Single lines represent covalent 
bonds. The parameter w is the distance between the incident 
hydrogen atom and the surface of the graphene. The double 
lines represent expanding covalent bonds, and the parameter 
v is the length of a side of the hexagonal hole. 



therefore, we substitute the relative momentum for po 
and can perform similar estimation to the preceding sub- 
sections. As a result, the absorption rate increases and 
the reflection rate decreases as the graphene temperature 
is raised. By comparison energy order, the penetration 
rate is insensitive to the graphene temperature. 



The graphene temperature dependence of reaction 
rates is significant at low incident energies in Fig. 2J As 
the graphene temperature is raised, the absorption rate 
increases and the reflection rate decreases for E\ <1 eV. 
The maximum temperature 2000 K, which corresponds 
to 0.26 eV as kinetic energy per a carbon atom, is com- 
parable to the energy height of the potential wall of the 
second repulsive force of 0.5 eV. If the nearest carbon 
atom has kinetic energy, the relative momentum between 
the hydrogen atom and the nearest carbon atom becomes 
larger than the initial momentum of the hydrogen atom 
Po in Eq. ([3]). In the case of high graphene temperature, 



V. SUMMARY 

By the CMD simulation with modified Brenner's 
REBO potential model, we demonstrated the chemical 
reaction between the single hydrogen atom and the single 
graphenc, which can be regarded as the elemental reac- 
tion between hydrogen and graphitic carbon materials. 
We observed the three processes, which are the absorp- 
tion, the reflection and the penetration reactions. The 
dominant reaction is replaced according to the incident 
energy for 0.1 eV < E\ < 100 eV. We discussed the 
characteristic interactions between the hydrogen atom 
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and the graphene with potential energy. The hydro- 
gen atom receives the repulsive force not only by nu- 
clei of carbon atoms but also by 7r-elcctrons over the 
surface of the graphene. These two kinds of repulsive 
force cause the two reflection mechanisms. When the 
hydrogen atom is absorbed, the graphene is transformed 
from flat sheet configuration to "overhang" configuration. 
By comparison between the typical time length of the 
overhang transformation and the time length during the 
hydrogen atom's stay, we can clarify the difference be- 
tween the absorption reaction and the reflection reaction 
for E\ > 0.5 eV. In the penetration reaction, the inci- 
dent hydrogen atom goes through the hexagonal hole of 
the graphene and the graphene expands the hexagonal 
hole, simultaneously. The expansion lowers the energy 
height of the potential wall by nuclei of the carbon atoms, 
which accounts for starting the penetration reaction at 
E\ = 13 cV in Fig. [3l In addition, we investigated the 
graphene temperature dependence of reaction rates. As 
the graphene temperature rises, the absorption rate in- 
creases and the reflection rate decreases for low incident 
energy. The cause of the graphene temperature depen- 
dence is that the kinetic energy of the nearest carbon 
atom is comparable to the energy height of the potential 
wall by 7r-electrons. 



segment which starts at the i-th atom and ends at the 
j-th atom and the line segment which starts at the i-th 
atom and ends at the fc-th atom, as follows: 



COS Of 



(A.2) 



where Xi is the position coordinate of the i-th atom and 
Tij is the distance between the i-th and the j-th atoms.. 
The dihedral angle 6*°-^ * s the angle between the triangle 
formed by the j-th, the i-th and the fc-th atoms and the 
triangle formed by the i-th, the j-th and the Z— th atoms. 

qDH 
7 kijl 



The cosine function of is given by 



/iDH ( X i x k) x ( x j x i) { x j Xi) X (xi Xj) 
COS tfr. 



kijl 



(A.3) 



The repulsive function VSi (r^ ) and the attractive 
function VjA.j (r ^ ) are defined by 

V m(^) = (l + ^) A M cxp (-a mnj ) , 

(A.4) 
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APPENDIX: MODIFIED BRENNER'S REBO 
POTENTIAL MODEL 

We note here the review of Brenner's reactive empiri- 
cal bond order (REBO) potential^ and our modification 
points. This potential model follows in the wake of Morse 
potential^ Abell potential 1 ^ and Tersoff potentialJ^^ 

The potential function U is defined by 



U 



E fe( r «)- 5 8(W 1 ^}>{« DH })^(r«) 



i,j>i 



(A.l) 

where r%j is the distance between the i-th and the j-th 
atoms. The bond angle 0^ ik is the angle between the line 



The square bracket such as [ij] means that each func- 
tion or each parameter depends only on the species of 
the i-th and the j-th atoms, for example V^, V^ H 
and V(? H (= V£ c ). The coefficients Q^, A {lj] , a [y ], 
B n [ij] and (3 n [ij] are given by Table U 

The cutoff function //^(fy-) determines effective 
ranges of the covalent bond between the i-th and the j- 
th atoms. Two atoms are bound with the covalent bond 
if the distance ry is shorter than -D^ 1 ] 11 - Two atoms are 
not bound with the covalent bond if the distance r,j is 
longer than D^ x . The cutoff function connects 
the above two states smoothly as 



r i 

i 

2 



1 + C0S(7T 



[ij] 



-D" 



r) 



(D™f < x < Dffi), 
(x > DgHp. 



(A.6) 

The constants £>™jj n and -D™j] X depend on the species of 
the two atoms (Table ITT|). The cutoff function /^(ry) 
distinguishes the presence of the covalent bond between 
the i-th and the j-th atoms. 

The potentials VSi and VjA, in Eq. (|A.1|) generate 
two-body force, because both are the function of the 
only distance r^. The multi-body force is used in- 
stead of the effect of an electron orbital. In this model, 
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TABLE I: The parameters 


for the repulsive function Vuji 


and the attractive function V£i. 


They depend on the species of the 


i-th and the j-th atoms. 








l l J I 




.Parameter 


V V 


MM 


V 'JTL Ol 1 1 V 


Q[ij] 


0.3134602960833 A 


0.370471487045 A 


0.340775728 A 


Ar-i 


10953.544162170 eV 


32 817355747 eV 


149 94098723 eV 


"Ml 


4.7465390606595 A -1 


3.536298648 A -1 


4.10254983 A -1 




12388.79197798 eV 


29.632593 eV 


32.3551866587 eV 


B 2 [ij] 


17.56740646509 eV 


eV 


eV 


B 3 [ij] 


30.71493208065 eV 


eV 


eV 




4.7204523127 A -1 


1.71589217 A -1 


1.43445805925 A -1 




1.4332132499 A -1 


oA- 1 


oA- 1 


03[ij) 


1.3826912506 A" 1 


OA" 1 


OA" 1 



TABLE II: The constants for the cutoff function f^Arij). TABLE IV: The parameters for the sixth order spline function 
They depend on the species of the i-th and the j— th atoms. 7c(cos6 lB 



"jikj 



gjf (A) (A) 



CC 1.7 2.0 

CH 1.3 1.8 

HH 1.1 1.7 



cos e% 


7c 


7c 


7c 


(3) 

7c 


cos(109.47°) 


0.09733 


0.400 


1.980 


-9.9563027 


1 


1.0 


0.78 


-11.3022275 





TABLE V: The parameters for the sixth order spline func- 
TABLE III: The parameters for the sixth order spline function tion G H (cos^ fe ). The parameters are determined under 



Gc(cOB0? tt ). 


cos 6% 


Gc 


G'c 


G'c 




-1 


-0.001 


0.10400 








-1/2 


0.05280 


0.170 


0.370 


-5.232 


cos(109.47°) 


0.09733 


0.400 


1.980 


41.6140 


1 


8.0 


0.23622 


-166.1360 





hj({r}, {6> B }, {6> DH }) in Eq. (|A~T|) gives multi-body force 
and is denned by 



MM,{0,{0) 



b-r({r},{0 B }) + b°r({r},{0 B }) 
U?f({r}) + b™({r},{9™}). 



(A.7) 



The first term ^ 



generates three-body force except 
the effect of 7r-electrons. The second term IT^ C in Eq. 
0A/7J) represents the influence of radical energetics and 
7r-bond conjugation^ The third term bf^ ({r} , {9 Dii }) 
in Eq. (|A.7|) derives four-body force in terms of dihedral 
angles. These functions are composed of the production 
of cutoff functions /^-j (Vjj ) ■ Five or more - body force 
are generated during chemical reaction. 

The function b^({r},{9 B }) in Eq. JX7]) is defined 

by 



b-ri{r},{0 B }) 



l + E fM(rr 3 )G t (cos9f ik )e' 

_ i 

+ P M (N^,Np j )] \ (A. 



cos 9 



jik 



Parameter 



Value 



G H (0) 
Gk(0) 
G H (0) 

4 3) (o) 

G«(0) 
Gk 5) (0) 
Gk 6) (0) 



19.06510 
1.08822 

-1.98677 
8.52604 

-6.13815 

-5.23587 
4.67318 



The function Gi in Eq. (|A.8j) depends on the species of 
the i-th atom. If cos9 B ik > cos(109.47°) and the i-th 

atom is carbon, G; is defined by 



Gi (cos 9f lk ) = [1 - Q c ( AfJ)] G c (cos 6%) 
+Q c (Ml) lc (co S ef ik ). 



(A.9) 



If cos 9 B ik <J cos(109.47°) and the i— th atom is carbon, 
Gi is defined by 

G^cos^JeGcIcos^). (A.10) 

And, if the i-th atom is hydrogen, Gi is defined by 

^(cos^JeeGh^os^,). 



(A.ll) 



Here Gc, 7c and Gh are the sixth order polynomial 
spline functions. Though the spline function G; needs 
seven coefficients, the only six coefficients are written in 
Brenner's paper £ We determine the seven coefficients in 
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TABLE VI: Parameters for the bicubic spline function 
P[ij](Nfj , N^). The parameters which are not denoted are 
zero. 





Value 


Pec (1,1) 


0.003026697473481 


Pcc(2,0) 


0.007860700254745 


Pec (3,0) 


0.016125364564267 


Pcc(l,2) 


0.003179530830731 


Pec (2,1) 


0.006326248241119 


Pch(1,0) 


0.2093367328250380 


Pch(2,0) 


-0.064449615432525 


Pch(3,0) 


-0.303927546346162 


Pch(0,1) 


0.01 


Pch(0,2) 


-0.1220421462782555 


Pch(1,1) 


-0.1251234006287090 


Pch(2,1) 


-0.298905245783 


Pch(0,3) 


-0.307584705066 


Pch(1,2) 


-0.3005291724067579 



table IIII1 IIVI and [Vj respectively. The function Q c and 
the coordination number M\ in Eq. (|A.9[) are defined by 



X*) 




(x < 3.2), 
cos (2?r (x - 3.2))] (3.2 < x < 3.7) , 
(x > 3.7), 

(A.12) 



E f[ik] 



ink)- 



(A.13) 



The constant Xujk] in Eq. (|A.8[) is a weight to mod- 
ulate a strength of three-body force, which depends on 
the species of the i-th, the j-th and the k-th atoms. 
In comparison with Brenner's former potential^ we set 
constants \ijk] as follows: 

Ahhh = 4.0, (A.14) 
Accc = Acch = Achc = Ahcc 

= Ahhc = Ahch = Achh = 0. (A. 15) 



The function P\ij] in Eq. (|A.8p is required in the case 
that molecules forms solid structure. The function P[ij\ 
is the bicubic spline function whose coefficients depend 
on the species of the i-th and the j-th atoms (Table 
IVl*)) . The parameters Njj and N$ are, respectively, the 
number of hydrogen atoms and the number of carbon 
atoms bound by the i-th atom as follows: 



hydrogen 

E /fofok). (A.16) 

carbon 

E /p*] (»■*)• ( A - 17 ) 

The second term II? C in Eq. (|A.7|) is defined by a 
tricubic spline function as 

n* c ({r}) ee F M (N^,N^N^), (A.18) 



TABLE VII: Parameters for the tricubic spline func- 
The parameters which are not denoted are 



tion 
zero. 



The function Fuj] satisfies the following rules: 
F [y] (JVi,JV a ,JVa) = F M (N a ,Ni,N a ), d Nl F M (N 1 ,N 2 ,N 3 ) = 
d Nl F [i:j] (N 2 ,Ni,N 3 ), F fo] (JVi,Ar 2 ,JV3) = F m (3, N 2 , N 3 ) if 
JVi > 3, and F M (N U N 2 ,N 3 ) = F lij] (N 1 ,N 2 ,5) if 7V 3 > 5, 
where (9jv ; = d/dNi. 



Variables 



Function 



Fcc(Ni,N 2 ,N 3 ) 



Fhh(N u N 2 ,N 9 ) 
Fca(Ni,Na,N s ) 



Ni 


N 2 


N 3 


Value 


1 


1 


1 


0.105000 


1 


1 


2 


-0.0041775 


1 


1 


3 to 5 


-0.0160856 


2 


2 


1 


0.09444957 


2 


2 


2 


0.04632351 


2 


2 


3 


0.03088234 


2 


2 


4 


0.01544117 


2 


2 


5 


0.0 





1 


1 


04338699 





1 


2 


0.0099172158 





2 


1 


0493976637 


o 


2 


2 


-0.011942669 


o 


3 


1 to 5 


-0.119798935 


1 


2 


1 


0.0096495698 


1 


2 


2 


0.030 


1 


2 


3 


—0 0200 


I 


2 


4 to 5 




1 


3 


2 to 5 


-0.124836752 


2 


3 


1 to 5 


-0.044709383 


2 


1 


1 


-0.052500 


2 


1 


3 to 5 


-0.054376 


2 


3 


1 


0.0 


2 


3 


2 to 5 


0.062418 


2 


2 


4 


-0.006618 


1 


1 


2 


-0.060543 


1 


2 


3 


-0.020044 


1 


1 


1 


0.249831916 





2 


3 to 5 


-0.0090477875161288110 


1 


3 


1 to 5 


-0.213 


1 


2 


1 to 5 


-0.25 


1 


1 


1 to 5 


-0.5 



where the variables are defined by 

N h = E /[*](^)> 



(A.19) 



carbon 



with 



N^p = i+ J2 /[ ife ](n fe )c N (iv^) 

carbon 

+ E fvqfaOOsWj, (A.20) 



1 



(x<2), 

C N (x)= { | [1 + cob(tt(x - 2))] (2<a;<3), 
(x > 3). 

(A.21) 

The second and the third terms of the right hand of Eq. 
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TABLE VIII: Parameters for the tricubic spline function 
Tec- The parameters which are not denoted are zero. The 
function Tec satisfies the following rule: Tcc(Ni, N2, N3) = 
T C c{N u N 2 ,5) if N 3 > 5. 





Variables 




Function 


Ni N 2 N 3 


Value 


Tac(Ni,N 2 ,Na) 


2 2 1 


-0.070280085 




2 2 5 


-0.00809675 



(|A.20[) are not squared. We note that they are squared 
in Brenner's original formulation.— By this modification, 
a numerical error becomes smaller than Brenner's forma- 
tion. Table I VIII shows the revised coefficients for Tjy] . 
The third term i£ H ({r}, {6> DH }) in Eq. (|A~7]) is defined 



by 



&S H (W> {# DH }) = T[ij] (Nij>Njii 



EE 



9 DH ) 



/p](nfc)/g 7 ](^0 



(A.22) 



where Tpyi is a tricubic spline function and has the same 
variables as .Fry] in Eq. (|A.18[) 



Tryi is also revised due 



The coefficients for 
to the modified A^" J (Table 
IVIIip . In the present simulation, the function Try] be- 
comes Tcc(2, 2, 5) for a perfect crystal graphene, and be- 
comes Tqc (2, 2, 3) or Tcc(2, 2, 4) when a hydrogen atom 
is absorbed. 
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